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Abstract - Networks with a prescribed power-law scaling in the spectrum of the graph Laplacian 
can be generated by evolutionary optimization. The Laplacian spectrum encodes the dynamical 
behavior of many important processes. Here, the networks are evolved to exhibit subdiffusive 
dynamics. Under the additional constraint of degree-regularity, the evolved networks display an 
abundance of symmetric motifs arranged into loops and long linear segments. Exploiting results 
from algebraic graph theory on symmetric networks, we find the underlying backbone structures 
and how they contribute to the spectrum. The resulting coarse-grained networks provide an 
intuitive view of how the anomalous diffusive properties can be realized in the evolved structures. 


Introduction. — Networks have become a principal 
tool for the modeling of complex systems in a broad range 
of scientific fields Bi In a dynamical network the topol¬ 
ogy describes the couplings between individual units of a 
dynamical process . The general underlying question of 
research on dynamical networks is how the network struc¬ 
ture shapes the global dynamical behavior. 

An important class of processes on networks are Lapla¬ 
cian dynamics in which the graph Laplacian is the opera¬ 
tor of the (linearized) time evolution. The Laplacian spec¬ 
trum and eigenvectors then describe the overall dynamical 
behavior. This class comprises many physical processes 
such as dynamics of Gaussian spring polymers [^, trans¬ 
port processes 6], r andom walks [^, and synchronization 
of oscillators [8 TO . 

Another dynamical aspect is the evolution of network 
structure. Many systems change their connectivity struc¬ 
ture in the course of time which affects the global dy¬ 
namical behavior. The time scales of these two processes 
are often well separated with fast dynamics and slowly re¬ 
sponding structural evolution. As an example, think of 
changes in neuronal activity in a brain on the one hand 
and the formation of synaptic connections on the other 
hand. In the opposite case, when dynamics and evolution 
happen on similar time scales one speaks of coevolution¬ 
ary or adaptive networks II . Since the functionality of a 


system is often closely associated with dynamical behav¬ 
ior evolutionary forces will drive an evolving dynamical 
network towards some optimized structure. 

The method of network evolution adopts this strategy 


to find networks with prescribed dynamical properties. 
Based on rules for mutation and selection together with 
a “fitness” function measuring the quality of a network 
structure the evolution explores the configuration space 
towards optimal structures. Evolutionary optimization 
of networks was applied successfully in various contexts 
such as modularity in changing environments 12 , Boolean 
switching dynamics 13 T8 , and synchronization of oscil¬ 


lators 19 -23 . It has also been shown that networks can 


be reconstructed from their Laplacian spectra by evolu¬ 
tionary optimization 


24,25 


Recently, network evolution was applied to generate net- 

Here, sub¬ 


works which display subdiffusive dynamics 26 


diffusion is specified by the average return probability of 
a random walker that decays more slowly than in normal 
diffusion. The return probability is directly related to the 
spectrum of the graph Laplacian via a Laplace transform. 
Thus, in contrast to other questions like the synchroniz- 
ability of oscillators where mainly the smallest eigenvalues 


contribute 10 , here the entire spectrum is important. The 


evolution algorithm optimizes the spectrum with respect 
to a prescribed target function, which in this case is a pure 
power law characterized by the spectral dimension dg. For 
subdiffusive dynamics dg <2. 

To elucidate how the spectral properties are encoded 
in the structure of the evolved networks, in the present 
work we apply the method of ref. 26 to regular net¬ 


works, i.e., networks in which each vertex has the same 
number of neighbors. Two typical network configurations 
evolved towards a target spectral dimension of dg = 1.4 
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Fig. 1: (Color online) Top: Typical fc-regular network config¬ 
urations after 10® evolution steps towards spectral dimension 
da = 1.4 starting from (a) a square lattice with A: = 4, = 361, 

M = 722 and (b) a honeycomb lattice with k = 3, N = 360, 
M = 540. Bottom: The corresponding s-quotients which are 
of sizes (c) N = 230, M = 294 and (d) N = 285, M = 337. In 
the enlarged segments, the colors indicate the orbits to which 
the respective vertices belong. 


are shown in fig. 0^) and (b). They exhibit an interest¬ 
ing, modular structure consisting of small symmetric mo¬ 
tifs arranged into fractal-like loops and linear segments of 
different lengths on larger scales. Given that subdiffusive 
transport is common in comb-like and self-similar struc¬ 
tures 27 , it is very tempting to attribute the subdiffusive 
character of the networks to this modular structure. We 
know, however, that primarily spectral and not structural 
properties determine the dynamical behavior [8 28 and 
that the spectrum of a network can, in general, not be 
constructed from the spectra of its subnetworks. Never¬ 
theless, in the present case the symmetry of the networks 
implies that this can be done in a systematic way. 


In the following we exploit concepts of algebraic graph 
theory to extract the large-scale structure of the evolved 
networks. This is achieved by identifying the symmetries 
associated with the small-scale motifs and constructing 
the corresponding quotient network. The quotient graph 
and in particular the corresponding simple graph represent 
a systematic coarse-graining of the original networks, and 


their spectra are considerably closer to the target function 
of the evolutionary algorithm. To set the stage for this 
investigation, we begin by providing some background on 
spectral graph theory and network symmetries. 

Spectra of dynamical networks. — A formal de¬ 
scription of a network with N vertices is given by the 
adjacency matrix A. It is e,n N x N matrix with ele¬ 
ments Aij = 1 if vertex i is connected with vertex j and 
Aij = 0 otherwise. For an undirected network A is sym¬ 
metric, Aij = Aji. In the case of a directed network one 
has to specify whether Aij describes the presence of the 
an edge from i to j or vice versa. Both conventions ex¬ 
ist in the literature. The graph Laplacian L is another 
matrix describing the structure of an undirected simple 
network (no self-loops, no multi-edges). It is defined as 
L = D — A where D is the diagonal matrix of vertex 
degrees Dij = kiSij. The degree ki = J2j ^ij of a ver¬ 
tex i is the number edges incident to the vertex. Note that 
this operator is sometimes called the algebraic Laplacian 
to distinguish it from similar operators such as different 
versions of the normalized Laplacian. 

Diffusion processes form a very important class of dy¬ 
namics. They describe, e.g., how a substance spreads in a 
medium or an opinion in a society. To be more specific, we 
consider a process where the flux of some substance along 
each edge in the network is proportional to the difference 
in the amount at the end vertices. The time evolution 
operator of this process on an undirected network is the 
graph Laplacian. An important global characteristic of 
such a process is the average probability Po{t) for a ran¬ 
dom walker to return to its starting point at time t. Po(t) 
is the Laplace transform of the Laplacia n eigen value den¬ 
sity p(A), Po{t) = dAexp(—At)/9(A) [^[^. Hence, if 

the integrated eigenvalue density /(A) = J dAp(A) scales 
as a power law, /(A) oc then Poit) decays as a power 

law, Po{t) oc The power-law exponent ds is the 

spectral dimension of the network. For regular lattices, it 
coincides with the Euclidean dimension 


29 


Since the quotient graph to be constructed below is di¬ 
rected, we need to generalize the graph Laplacian to di¬ 
rected networks. For this purpose we consider a diffusion 
process on a directed network, possibly containing multi¬ 
edges and self-loops. Let Xi describe the amount of the dif¬ 
fusing quantity on vertex i. The flux of this quantity along 
each outgoing edge should be proportional to Xi. Then Xi 
decreases at rate cXi for each outgoing edge i ^ j, and 
increases at rate cxj for each incoming edge j ^ i where 
c is a diffusion coefficient. Following the convention that 
Aij is given by the number of directed edges from i to j 
this yields the set of coupled differential equations 

Xi = AjiCXj — AijCXi = AjiXj — cxik°'^^ 


)^ 3 — 


°'^^x ■ 

'31 -^3 


( 1 ) 


where fc™* = J2j ^ij denotes the out-degree of vertex i 
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and 

- A,, (2) 

are the elements of the out-degree Laplacian. So, when 
considering diffusion processes on a network the out- 
degree Laplacian is the correct generalization of the graph 
Laplacian as time evolution operator to directed networks. 

Symmetries in networks. — In a symmetric net¬ 
work, the spectrum and, correspondingly, the set of eigen¬ 
vectors can be split into two classes. Some of the eigen- 
pairs stem from the symmetric motifs and are inherited 
by the network as a whole while the remaining ones are 
the eigenpairs of the underlying backbone structure, the 
quotient network. In order to make this statement more 
precise and formally tractable we first summarize some 
basic notions on the mathematical description of network 
symmetry. More details, exact definitions, and proofs can 
be found in ref. [M|[3T] . 

The symmetry of a network G is described by its auto¬ 
morphism group Aut(G'). An automorphism is a permuta¬ 
tion of vertices which does not change the network struc¬ 
ture (preserves the adjacency). A network is called sym¬ 
metric if it has a non-trivial automorphism group. By a 
direct product decomposition Aut(G') = H 1 XH 2 X .. .xHk 
the automorphism group can be split into geometric fac¬ 
tors Hi which act independently on the network. A sym¬ 
metric motif A4 h of a network is the induced subgraph 
on the support of a geometric factor H. Intuitively speak¬ 
ing, a symmetric motif Af// is a minimal subgraph whose 
vertices are moved by H. 

The action of a network’s automorphism group Aut(G') 
partitions the vertex set V into disjoint equivalence classes 
called orbits. The orbit of a vertex r; S F is the set of 
vertices A.(v) C 1^ to which v maps under the action of 
Aut(G), 

= { 5 ^^ G ^ 1 5 € Aut(G)} . (3) 

In the enlarged parts of fig. [^a) and (b) examples of how 
the vertices in a symmetric motif are grouped into orbits 
are indicated by their colors. In general, it may, however, 
be necessary to swap a larger fraction of a network in order 
to map two vertices in the same orbit onto each other 
without altering the adjacency. Nevertheless all vertices 
in the same orbit are structurally equivalent (perform the 
same function in the network). 

The quotient graph is a network coarse-graining based 
on this structural equivalence. The set of vertices of 
the quotient Q = G/A is given by the orbits A = 
{Ai, A 2 ,..., Aj.} of the network G. Note that due to the 
structural equivalence the number of neighbors in Ay of a 
vertex v G A^, called g^y, does not depend on the choice 
of V but only on x and y. A partition with this property 
is called equitable which has important consequences for 
the relation between the spectra of G and Q (see below). 
The quotient graph Q of a network G has vertex set A 
and adjacency matrix {qxy}x,y=i,...,r- It is by definition di¬ 
rected and generally contains multi-edges and self-loops. 


If one is working with simple networks this may be incon¬ 
venient. As a further simplification the s-quotient Qs is 
defined as the underlying simple graph of the quotient. It 
already captures essential properties such as the size and 
adjacency of the quotient. Fig. [^c) and (d) show the s- 
quotients of the two networks above. The colored vertices 
in the enlarged segments indicate the corresponding orbits 
in the parent networks. 

A key result from algebraic graph theory relates the 
adjacency spectra and eigenvectors of a network G and its 
quotient graph Q (or, in general, any quotient G/tt from 
an equitable partition tt). All eigenvalues of the quotient 
are also eigenvalues of the parent network. Moreover, to 
each eigenpair (A, v = (ui,... ,Vm)) of Q there exists an 
eigenpair (A, v = (iii,..., v„)) of G with Vi = Vx for all i G 
Ax: (eigenvectors are constant Vx on each orbit A^,). These 
eigenvectors are called lifted from the eigenvectors of the 
quotient. The remaining eigenpairs of G are constructed 
from the redundant eigenpairs of the symmetric motifs of 
G: If (A, V = (ui,..., Vk)) is an eigenpair of the symmetric 
motif Ai —isolated from the rest of the network—which is 
redundant (i.e. ~ ^ for each orbit A G M.) then 

(A,v) is an eigenpair of G with iii = Vi for a\\ i G Ai 
and Di = 0 for i ^ Af (eigenvectors are localized on the 
symmetric motifs). 

To summarize, the adjacency eigenpairs of a symmetric 
network G fall into two classes. (1) The eigenpairs of the 
quotient of G are also eigenpairs of G with constant val¬ 
ues of the eigenvectors on the orbits. (2) The redundant 
eigenpairs of the symmetric motifs of G are also eigenpairs 
of G with eigenvectors localized on the respective motif. 

All these statements for the adjacency spectrum and 
eigenvectors are elaborated and proven in ref. 30 . In the 


appendix we show that they also hold for the Laplacian 
spectrum and eigenvectors. 

For practical computations of automorphism groups ef¬ 
ficient and easily applicable algorithms are available such 


as the nauty program 32 which was used in this study 


to find orbits and construct network quotients. 

Application to regular evolved networks. — 

ref. 


In 


26 the method of network evolution was applied to 


construct networks with a given non-trivial spectral di¬ 
mension ds- The algorithm searches the configuration 
space of connected networks with a given size for network 
structures which resemble the power law in the Laplacian 
spectrum as closely as possible. The introduced distance 
measure V (called A in ref. [^) is defined as the inte¬ 
gral over the squared difference between integrated eigen¬ 
value densities. In contrast to other proposed spectral dis¬ 


tances 33 it can be applied to compare the spectrum with 


a general function as well as with other spectra. However, 
since T) inherently depends on the number of vertices, net¬ 
works of different sizes cannot be compared directly. The 
most simple rules for mutation and selection were applied, 
namely a random movement of a single edge and the ac¬ 
ceptance of a proposed move if T) decreases and the whole 
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Fig. 2: (Color online) Redundant Laplacian eigenvectors of the 
most prominent motifs in 4-regular (left) and 3-regular (right) 
evolved networks as depicted in fig. [^a) and (b). The sym¬ 
metric vertices are colored red. The corresponding eigenvalues 
are {5,5} (left) and {4} (right). 



network remains connected. Starting from regular lattices 
and random graphs networks were successfully evolved un¬ 
der the constraints of a fixed number of vertices and a fixed 
average degree. 

An even stricter condition is to fix the individual de¬ 
gree of each vertex, as for example in a fc-regular network 
where all vertices have the same degree k. Is it—under 
this additional constraint—still possible to construct net¬ 
works with a given non-trivial value of the spectral dimen¬ 
sion ds? To answer this question we applied the network 
evolution algorithm starting from square lattices {k = 4) 
and honeycomb lattices (fc = 3) and corresponding 4- and 
3-regular random networks. As in the preceding study [26] 
it turns out that the evolution does not depend on the two 
choices (lattice or random network) of initial conditions. 
The conservation of the degree sequence (or the individual 
degree of each node) in the course of the evolution can be 
realized by an “edge-crossing” update. In each evolution 
step two edges Vi-V 2 and V 3 -V 4 are randomly chosen and 
“crossed” to V 4 -V 4 and U 3 -U 2 . When the edges are drawn 
one has to make sure that the four vertices z;i,U 2 ,U 3 ,U 4 
are all distinct and that the edges V 1 -V 4 and V 3 -V 2 do not 
exist before the update. 

The evolutions were run for 10® time steps towards a 
target value of the spectral dimension of ds = 1.4. Typi¬ 
cal network configurations at the end of the evolutions of 
3- and 4-regular networks are shown in fig. Sa) and (b). 
These networks appear to have similar large-scale struc¬ 
tures of long loops built up by very frequently appearing 
symmetric motifs on small scales. The most abundant 
of these motifs are shown in fig. They consist of a 
complete subnetwork oi k — 1 vertices, each of which is 
joined to two boundary vertices connecting to the rest of 
the network. The vertices are labeled by their redundant 
Laplacian eigenvectors and the corresponding eigenvalues 
are {5, 5} in the 4-regular and {4} in the 3-regular case. 
These are also eigenvalues of the network as a whole, which 
explains the high degeneracy giving rise to the large step in 
the integrated eigenvalue density in fig. The motifs are 
assembled to form large-scale structures with long loops 
which appear to be quite similar in both cases. The large- 
scale structures of these networks become more evident in 


Fig. 3: (Color online) Logarithmically integrated Laplacian 
eigenvalue densities of the 4-regular evolved network from 
a) and its quotient, s-quotient (fig. [^c)) and an evolved 
network of the same size (number of vertices and number of 
edges) as the s-quotient. The eigenvalue density is plotted as 
a function of In A, where A = A/Amax < 1 is the eigenvalue 
normalized by the largest eigenvalue of the network 26 . The 


black dotted line indicates the target integrated density with 
d, = 1.4. 


their s-quotients depicted in fig. Be) and (d). One can 
clearly see the arrangement of loops of different lengths. 
Seen as a coarse-graining transformation, the s-quotient 
on average reduces the network sizes from = 361 ver¬ 
tices and M = 722 edges to N Pi 229(6) and M « 302(17) 
for 4-regular parent networks and from N = 360, M = 540 
to N Pi 281(6), M Pi 330(9) for 3-regular parent networks 
(standard deviations in parentheses). 

How to assess the quality of the s-quotients with re¬ 
spect to the optimization goal of the evolution? The 
ultimate goal of the evolution algorithm is to construct 
networks with a Laplacian spectrum which resembles the 
power law /(A) oc Figure displays the spectra— 

as logarithmically integrated densities /(In A) (see eq. (5) 
of ref. for the exact definition)—for an exemplary 4- 
regular evolved network, its quotient, s-quotient, and an 
evolved network of the same size as the s-quotient (“mixed 
and evolved s-quotient”). The latter class of networks was 
analyzed in order to get an idea of how the quotients and 
s-quotients compare to unconstrained evolved networks of 
the same size. They are constructed by optimizing net¬ 
works of the same size (same numbers of vertices and 
edges) by means of the unrestricted version of the evo¬ 
lution algorithm (conserving these numbers and the con¬ 
nectedness only). Technically, since the s-quotients are 
very sparse networks, instead of generating random net¬ 
works, the edges of the s-quotients are randomized (keep¬ 
ing the whole network connected) and then the evolution 
is run. 

As predicted by theory, the spectrum of the quotient re¬ 
sembles the spectrum of its parent up to the highly degen¬ 
erate redundant eigenvalues of the symmetric motifs. By 
removing the corresponding steps in the integrated eigen- 
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Fig. 4: (a) Typical 3-regular network configuration after 10® 
evolution steps towards ds = 1.1 starting from a honeycomb 
lattice with N = 360, M = 540 and (b) corresponding s- 
quotient, N = 260, M = 263. 


value density the fit to the target spectrum is improved 
in the regime of large eigenvalues, whereas the region of 
small A is unaffected. In contrast, there is no rigorous re¬ 
lation between the spectra of the parent network and the 
s-quotient. Therefore it was not to be expected that the 
spectrum of the s-quotient resembles the power law even 
more closely, as seen in the figure. In terms of the aver¬ 
age distance measure T) estimated from 100 realizations of 
the evolution, the quotients have a value of I? = 0.09(1) 
which is reduced by a factor of three in the s-quotients 
(T) = 0.03(2)). Nevertheless, the unconstrained evolution 
of networks of the same size succeeds in finding networks 
with a spectrum even closer to the desired power law, 
reaching an average value oi T) = 0.0045(1). However, 
similar to the networks obtained in ref. [^, the networks 
resulting from the unconstrained evolution algorithm can¬ 
not be easily spread out and visualized in a plane. 

What happens for a more extreme evolution tar¬ 
get? Since the evolution runs were started from two- 
dimensional lattices it should be more difficult to find net¬ 
works with a spectral dimension closer to 1. Fig. shows 
a typical network configuration from an evolution towards 
ds = 1.1 together with the corresponding s-quotient. In 
this case, the basic motifs are arranged to form even longer 
linear chains instead of loops resulting in an almost one¬ 
dimensional s-quotient. 

An important question is why the evolved networks as- 



Fig. 5: (Color online) Histogram of linear segment lengths s in 
s-quotients from 4-regular and 3-regular evolved networks for 
different target spectral dimensions dg. Shown are cumulative 
histograms H{s) = f°° h{t) dt of the absolute frequencies h(s) 
for 100 realizations of each. 


sume the observed shape of symmetric motifs assembled 
in loops and linear chains (see fig. [^a) and (b) andjijja)). 
A possible explanation is that they have too many edges 
to form structures with a low spectral dimension. The 
formation of the symmetric motifs makes it possible to 
“hide” the excess edges in the motifs and form large-scale 
structures with the desired scaling in the small Laplacian 
eigenvalues. Additionally, the motifs may effectively act 
as particle traps and thus slow down the diffusion process. 
The formation of loops and linear chains on larger scales is 
brought out more clearly in the coarse-grained s-quotients 
in fig. [^c) and (d) and |^b) where the symmetric motifs 
are removed. As a first approach towards quantifying the 
loop structure we analyzed the lengths of linear segments 
(connected maximal subgraphs of nodes with a maximum 
degree two) in the s-quotients in fig. We observe a very 
broad distribution in this quantity. However, the networks 
appear to be too small for a systematic scaling analysis. 


Conclusions. — We presented and analyzed subdiffu- 
sive networks found by evolutionary optimization. Subdif¬ 
fusion is known to be observed in various network struc¬ 
tures such as regular fractals, percolation clusters or combs 
with power-law distributed teeth lengths 27 . In the lat¬ 
ter, the delay is caused by the time a random walker 
spends on the teeth of the comb. In the networks pre¬ 
sented here, we observe backbone structures with loops 
and dangling ends of different length scales which seem to 
be the main source of subdiffusive behavior. The networks 
are, however, too small for a quantitative analysis of the 
distribution of these lengths. A more in-depth study on 
this question would be an interesting continuation of this 
work. 


* * * 

We thank Thilo Gross for sharing his knowledge on 


P-5 

















Steffen Karalus Joachim Krug 


network symmetries, Markus Porto for his support in 
early stages of this work, and Sungmin Hwang for valu¬ 
able comments and discussions. We gratefully acknowl¬ 
edge partial funding by the Studienstiftung des deutschen 
Volkes and the Bonn-Cologne Graduate School of Physics 
and Astronomy. 


Appendix. — We show that the relations between ad¬ 
jacency spectra and eigenvectors of a network G on the 
one hand and its quotient Q and symmetric motifs on 
the other hand also hold for the Laplacian spectrum and 
eigenvectors. The reasoning is the same as in Appendix A 
of ref. 30 . First, recall that the characteristic matrix P 
of the equitable partition tt = {Ci,... ,Cr} is defined by 
Pix = 1 if vertex i G Cx and Pix = 0 otherwise. It re¬ 
mains to prove that the space W spanned by the column 
vectors of P is also L-invariant. For this, W is L-invariant 
{Lu € W for all M G W) if and only if there exists a 
matrix L such that 


LP = PL. (4) 

L is then the graph Laplacian of the quotient network 
Q = G/tt. Let B {Bxy = qxy for the orbit partition above) 
be the adjacency matrix of the quotient. Define L = D—B 
with Dxy = kxSxy and kx = Pxy Since AP = PB it 
remains to show that DP = PD. The left hand side reads 

{DPfx — ^ ^ DijPjx = ^ ^ kiSijPjx — kiPix (5) 

j 3 

which is the degree of vertex iiii £ Gx and zero otherwise. 
The right hand side reads 

{.PP')ix ~ ^ ^ PiyPyx ~ ^ ^ PiykySyx 
y V 

— Pixkx — Pix ^ ^ Pxy • ( 6 ) 

V 

Since Bxy is the number of neighbors in Gy of a vertex 
from Cx, J2y Pxy gives the total number of neighbors of a 
vertex in Gx and eq. Q evaluates to the degree of vertex i 
if it belongs to Cx and zero otherwise. This is the same 
as eq. ([^ which completes the proof. 

Note that the definition of the Laplacian L of the 
quotient graph is consistent with the out-degree Lapla¬ 
cian defined in eq. ([^. In particular, the definition 
of Qxy corresponds to the convention that A^j describes 
directed edges from i to j. 
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